home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 8: LINUX Games / Linux Cubed Series 8 - LINUX Games.iso / games / video / fly8111-.000 / fly8111- / fly8 / oxplane.c < prev    next >
C/C++ Source or Header  |  1979-12-31  |  17KB  |  600 lines
  1. /* --------------------------------- oxplane.c ------------------------------ */
  2.  
  3. /* This is part of the flight simulator 'fly8'.
  4.  * Author: Eyal Lebedinsky (eyal@ise.canberra.edu.au).
  5. */
  6.  
  7. /* Dynamics of the eXpermental plane.
  8.  * This model will eventualy do a proper aerodynamic simulation. It will NOT
  9.  * go as far as non-linear 6DOF but will properly implement the linearized
  10.  * model. This program is being modified regularly as I learn more about
  11.  * the subject and gather more data. If I get a stable version then it will
  12.  * be separated into it's own module so as not to constantly have a half
  13.  * working program. But not yet.
  14.  *
  15.  * ToDo:
  16.  * 1 document units/scale used for all variables.
  17.  * 2 select aero formulas and break into classes:
  18.  *   - input parameters (from .prm)
  19.  *   - calculated fixed parameters
  20.  *   - slow changing parameters
  21.  *   - dynamically calculated variables.
  22.  * 3 use proper trans- and super-sonic calcs for Cd.
  23.  * 4 add wing sweep to calcs?
  24. */
  25.  
  26. #include "plane.h"
  27.  
  28.  
  29. /* Formulas used:
  30.  *
  31.  * b    wing span
  32.  * S    wing area
  33.  * ClRate0 slope of Cl for foil
  34.  * e    Oswald's (or span) efficiency factor
  35.  * aoa0    angle of attack where Cl is zero
  36.  * feff    degrees aoa effective for 1 degree flaps
  37.  *
  38.  * AR    wing aspect ratio
  39.  * k_AR    1/(pi*e*AR)
  40.  * ClRate  slope of Cl for wing
  41.  * e1    drag span-efficiency-factor. We use 'e' for it.
  42.  *
  43.  * V    airspeed
  44.  * rho    air density
  45.  * aoa    angle of attack (geometric)
  46.  * aoaf    angle of attack effect due to flaps
  47.  * aoai    induced angle of attack
  48.  *    = Cl / (pi*e1*AR), but folded into ClRate
  49.  * aeff    effective angle of attack
  50.  *  ... all of these angles in degrees!
  51.  * q    dynamic pressure
  52.  * Cl    lift coefficient
  53.  * Cd    drag coefficient
  54.  * Cdi    induced Cd
  55.  * Cdp    profile Cd
  56.  *
  57.  * AR = b^2 / S
  58.  * aoai = Cl / (pi*e1*AR)
  59.  * ClRate = ClRate0 / (1+57.3*ClRate0/(pi*e1*AR))
  60.  *
  61.  * aoaf = flaps*feff
  62.  * aeff = aoa + aoaf
  63.  * Cl = (aeff - aoa0) * ClRate
  64.  *
  65.  *
  66.  *
  67.  *
  68.  * Cdi = Cl^2 /(pi*e*AR)
  69.  * Cdp = Cdp_characteristic + Cd_flaps + Cd_gear + Cd_bombs + Cd_speedbrakes
  70.  * Cd = Cdi + Cdp
  71.  * q = rho * V^2 / 2
  72.  * lift is perpendicular to airflow:
  73.  *     lift = Cl * q * S
  74.  * drag is parallel to airflow:
  75.  *     drag = Cd * q * S
  76.  *
  77.  * We need these factors:
  78.  * Ixx    pitching moment of inertia
  79.  * Iyy    rolling moment of inertia
  80.  * Izz    yawing moment of inertia
  81.  * we calculate the moments as:
  82.  * Mx = lift*Lc + drag*Dc + thrust*Tc + (stablilizer+elevators)*Sc
  83.  * Lc, Dc, Tc, Sc: effective moments levers.
  84.  * My,Mz = f(many things...)
  85. */
  86.  
  87. extern void FAR
  88. dynamics_xplane (OBJECT *p, int action)
  89. {
  90.     POINTER    *ptr;
  91.     int    onground, taxiing;
  92.     int    t, vmag, weight, thrust, push, drag, lift, k_ar;
  93.     int    rho, sos, Cd, Cl, ClTail, Cdi, sa, ca, sb, cb;
  94.     int    aFlaps, ClMax, ClMin, ClStall, Clrate, stall;
  95.     int    rVS, S, b2, b2VV, c2, c2VV, vmagV;
  96.     int    Crudder, Cailerons, Celevators, Cb, Cdx, Cdy, Cdz;
  97.     int    Fx, Fy, Fz, Mx, My, Mz, Ixx, Iyy, Izz;
  98.     long    tt;
  99.     ANGLE    alpha, beta, a;
  100.     VECT    AA;
  101.     AVECT    da;
  102.  
  103.     if (action)
  104.         return;
  105.  
  106.     if (dynamics_input (p))
  107.         return;
  108.  
  109.     ptr = p->pointer;
  110.  
  111.     onground = EX->flags & PF_ONGROUND;
  112.  
  113. /* keep the load within the designate range by limiting the elevators.
  114. */
  115.     if ((t = EX->Gforce - (EP->max_lift-EP->max_lift/4)) > 0) {
  116.         if (t > EP->max_lift>>1)
  117.             EX->elevators = 0;
  118.         else
  119.             EX->elevators -= muldiv (EX->elevators, t,
  120.                             EP->max_lift>>1);
  121.     } else if ((t = EX->Gforce - (EP->min_lift-(EP->min_lift>>2))) < 0) {
  122.         if (t < EP->min_lift>>1)
  123.             EX->elevators = 0;
  124.         else
  125.             EX->elevators -= muldiv (EX->elevators, t,
  126.                             EP->min_lift>>1);
  127.     }
  128.     EX->flags &= ~PF_GLIMIT;
  129.  
  130.     AA[X] = -fmul (GACC, p->T[X][Z]);    /* forces: gravity */
  131.     AA[Y] = -fmul (GACC, p->T[Y][Z]);
  132.     AA[Z] = -fmul (GACC, p->T[Z][Z]);
  133.  
  134.     EX->Gforce = -AA[Z];            /* pilot's gravity */
  135.  
  136.     da[X] = da[Y] = da[Z] = 0;        /* rotations */
  137.  
  138.     if (onground ? check_takeoff (p) : check_land (p)) {
  139.         p->flags |= F_HIT;
  140.         return;
  141.     }
  142.  
  143.     taxiing = onground && p->speed < 40*VONE;    /* NWS */
  144.  
  145.     airdata (p->R[Z], 0, 0, &rho, &sos);
  146. CCshow (p, 3, "rho", (long)fmul(rho, 1000));
  147. CCshow (p, 0, "sos", (long)DV(sos));
  148.  
  149. /* Get engine thrust.
  150. */
  151.     f16engine (p, sos);
  152.  
  153. /* Automatically refuel when stationery on ground.
  154. */
  155.     if (onground && 0 == p->speed && 0 == EX->power)
  156.         supply (p, 0);
  157.  
  158.     thrust = fmul (EX->thrust*2, FCON(0.453*9.8/VONE*5));    /* N*VONE */
  159.     tt = EP->weight + EX->fuel/100;                /* lb */
  160.     weight = fmul ((int)DV(tt), FCON (0.453));        /* Kg*VONE */
  161.     if ((t = EX->stores[WE_M61-1]) > 0)
  162.         weight += (int)DV(t * st.bodies[O_M61]->shape->weight/1000);
  163.     if ((t = EX->stores[WE_MK82-1]) > 0)
  164.         weight += (int)DV(t * st.bodies[O_MK82]->shape->weight/1000);
  165.     push = muldiv (thrust, VONE, weight);
  166. CCshow (p, 0, "thrust", (long)thrust);
  167. CCshow (p, 0, "weight", (long)weight);
  168. CCshow (p, 0, "thrA", (long)push);
  169.  
  170.     if (taxiing) {
  171.  
  172.         EX->flags &= ~(PF_GLIMIT|PF_STALL);
  173. /* Taxiing, ignore aerodynamics altogether. Only NWS and wheels friction.
  174. */
  175. CCshow (p, 0, "vold", (long)p->speed);
  176.         p->speed += TADJ(push);            /* v */
  177. CCshow (p, 0, "vnew", (long)p->speed);
  178.         drag = EP->wheel_mu
  179.             + muldiv (EP->brake_mu-EP->wheel_mu, EX->brake, 100);
  180.         drag = TADJ (fmul (drag, GACC));    /* assume no lift */
  181. CCshow (p, 0, "drag", (long)drag);
  182.         if (p->speed > drag)
  183.             p->speed -= drag;
  184.         else if (p->speed < -drag)
  185.             p->speed += drag;
  186.         else
  187.             p->speed = 0;            /* stop! */
  188. CCshow (p, 0, "v", (long)p->speed);
  189.  
  190.         if (p->a[X] > EP->gpitch)
  191.             p->a[X] -= TADJ(VD90);        /* slam down */
  192.         if (p->a[X] < EP->gpitch)
  193.             p->a[X] = EP->gpitch;
  194.  
  195.         if (p->a[Y])
  196.             p->a[Y] -= TADJ(p->a[Y]);    /* level wings */
  197.  
  198.         EX->v[X] = 0;
  199.         EX->v[Y] = fmul (p->speed,  p->cosx);
  200.         EX->v[Z] = fmul (p->speed, -p->sinx);
  201.  
  202.         p->da[X] = p->da[Y] = 0;
  203.         t = muldiv (EP->MaxRudder, EX->rudder, RAD2ANG(100));
  204.         p->da[Z] = muldiv (p->speed, t, VONE); /* NWS */
  205.         da[X] = da[Y] = da[Z] = 0;
  206.     } else {
  207.  
  208. /* Flying, although may still be on ground. Must do full arerodynamics but
  209.  * allow for ground contact.
  210.  *
  211.  * Ideal Cl rate is 2*pi, and for a given aspect ratio A it is 2*pi*A/(A+2).
  212.  * We represent angles as 2.0 angle = pi radians, so the formula is:
  213.  * pi*pi*A/(A+2). We prefer to keep the inverse which is: (1+2/A)/(pi*pi).
  214.  * And to avoid fraction overflow inside the parens: 2*(0.5+1/A)/(pi*pi). We
  215.  * keep 1/A in k_ar [which later is 1/(pi*A)].
  216. */
  217.         k_ar = fdiv (EP->wing_area,
  218.                 muldiv (EP->wing_span, EP->wing_span, VONE));
  219. CCshow (p, 3, "1/AR", (long)fmul(k_ar, 1000));
  220.         Clrate = fmul (FCON(2.0/(C_PI*C_PI)), FCON(0.5) + k_ar);
  221. CCshow (p, 2, "Clrate", (long)ANG2DEG00 (Clrate));
  222.         k_ar = fmul (k_ar, FCON(1.0/C_PI));        /* 1/(pi*AR) */
  223. CCshow (p, 3, "k/AR", (long)fmul(k_ar, 1000));
  224.  
  225.         alpha = EX->v[Y] ? -ATAN (EX->v[Z], EX->v[Y]) : 0;
  226.         alpha += EP->Aoffset;
  227.         beta = p->speed ? ASIN (fdiv (EX->v[X], p->speed)) : 0;
  228. CCshow (p, 2, "alpha", (long)ANG2DEG00 (alpha));
  229. CCshow (p, 2, "beta", (long)ANG2DEG00 (beta));
  230.         sa = SIN (alpha);
  231.         ca = COS (alpha);
  232.         sb = SIN (beta);
  233.         cb = COS (beta);
  234.         EX->aoa = alpha;
  235.  
  236.         a = alpha - EP->Cl0;
  237.         ClTail = muldiv (EP->Tvol, a+EP->Toffset, Clrate);
  238. CCshow (p, 3, "ClTail", (long)fmul(ClTail, 1000));
  239.  
  240.         aFlaps = fmul (EP->FEff, muldiv (EP->MaxFlaps, EX->flaps, 100));
  241.         a += aFlaps;
  242. CCshow (p, 2, "aeff", (long)ANG2DEG00 (a));
  243.  
  244. /* 'CLf' is to avoid overflow (fractions are limited to the range [-2...+2)).
  245. */
  246. #define    CLf    4
  247.         if (iabs(a>>1) < (Uint)(CLf*Clrate))
  248.             Cl = muldiv (FONE/CLf, a, Clrate);
  249.         else if (a > 0)
  250.             Cl =  FCON(1.99);
  251.         else
  252.             Cl = -FCON(1.99);
  253.  
  254.         t = muldiv (FONE/CLf, aFlaps, Clrate);
  255.         ClMax = EP->maxCl/CLf*10 + t;
  256.         ClMin = EP->minCl/CLf*10 + t;
  257.         EX->flags &= ~PF_STALL;
  258.         stall = 0;
  259.         ClStall = FONE;
  260.         if (Cl > ClMax) {
  261.             EX->flags |= PF_STALL;
  262.             if (!(EX->flags & PF_NOSTALL)) {
  263.                 stall = 1;
  264.                 Cl = Cl/4 - (Cl-ClMax);
  265.                 if (Cl < 0)
  266.                     Cl = 0;
  267.                 else
  268.                     Cl *= 4;
  269.             } else {
  270.                 Cl = ClMax;
  271.                 ClStall = iabs (ca);
  272.             }
  273.         } else if (Cl < ClMin) {
  274.             EX->flags |= PF_STALL;
  275.             if (!(EX->flags & PF_NOSTALL)) {
  276.                 stall = 1;
  277.                 Cl = Cl/4 + (ClMin-Cl);
  278.                 if (Cl > 0)
  279.                     Cl = 0;
  280.                 else
  281.                     Cl *= 4;
  282.             } else {
  283.                 Cl = ClMin;
  284.                 ClStall = iabs (ca);
  285.             }
  286.         }
  287.  
  288.         Cdi = fmul (Cl, k_ar*CLf);
  289.         Cdi = fmul (Cdi, Cl)*CLf;
  290.         Cdi = muldiv (Cdi, 100, EP->efficiency_factor);
  291.  
  292. /* The groung effect formula is extracted from the graph in Smiths 'The
  293.  * Illustrated Guide To Aerodynamic' end of chapter 3.
  294. */
  295.         t = EP->wing_span;
  296.         if (p->R[Z] < (long)t) {        /* ground effect */
  297.             t = fdiv ((int)p->R[Z], t);    /* h/b */
  298. CCshow (p, 3, "h/b", (long)fmul(t, 1000));
  299.             if (t < FCON(0.1))
  300.                 t = 5 * t;
  301.             else
  302.                 t = FCON(1.06) - fdiv (FCON(0.07), t);
  303. CCshow (p, 3, "gef", (long)fmul(t, 1000));
  304.             Cdi = fmul (Cdi, t);
  305.         }
  306. CCshow (p, 3, "Cl", (long)fmul(Cl, CLf*1000));
  307. CCshow (p, 3, "Cdi", (long)fmul(Cdi, 1000));
  308.  
  309.         Cd = EP->Cdp0;
  310. CCshow (p, 3, "Cdp0", (long)fmul(Cd, 1000));
  311.         Cd += Cdi;
  312. CCshow (p, 3, "+Cdi", (long)fmul(Cd, 1000));
  313.         if (EX->airbrake) {
  314.             Cd += muldiv (EP->Cds, EX->airbrake, 100);
  315. CCshow (p, 3, "+Brks", (long)fmul(Cd, 1000));
  316.         }
  317.         if (EX->equip & EQ_GEAR) {
  318.             Cd += EP->Cdg;
  319. CCshow (p, 3, "+Gear", (long)fmul(Cd, 1000));
  320.         }
  321.         if (EX->stores[WE_MK82-1] > 0) {
  322.             Cd += EX->stores[WE_MK82-1] * EP->CdMK82;
  323. CCshow (p, 3, "+MK82", (long)fmul(Cd, 1000));
  324.         }
  325.  
  326. /* We now calculate the new forces, based on the state of the plane.
  327.  * The given state is:
  328.  *  Cdx        Rotation around x (nose up, pitch)
  329.  *  Cdy        Rotation around y (right wing down, roll)
  330.  *  Cdz        Rotation around z (nose right, yaw)
  331.  *  Crudder    rudder right-turn angle
  332.  *  Cailerons    ailerons roll-right angle
  333.  *  Celevators    elevators pull-up angle
  334.  *  alpha    nose above wind angle
  335.  *  Cb        nose left-of-wind angle
  336.  * The forces are:
  337.  *  Fx        right
  338.  *  Fy        forward
  339.  *  Fz        up
  340.  * The moments are:
  341.  *  Mx        around x (nose up, pitch)
  342.  *  My        around y (right wing down, roll)
  343.  *  Mz        around z (nose right, yaw)
  344.  *
  345.  * From each force we get acceleration [a = F/m] and then velocity [v += a*t]
  346.  * and position [x += (v + a*t/2)*t].
  347.  *
  348.  * From each moment we get angular acceleration [aa = M/I], then we get
  349.  * rotation rate [da += aa*t] and angle [a += (da + aa*t/2)*t].
  350. */
  351.         vmagV = p->speed;
  352.         vmag = DV(vmagV);
  353.  
  354. /* 0.093 is for converting 'ft^2' to 'm^2'
  355. */
  356.         S = DV(EP->wing_area);
  357.         rVS = muldiv (fmul (fmul (rho, FCON(1.225)), vmag), S, VONE);
  358. CCshow (p, 0, "rVS", (long)rVS);
  359.         b2 = EP->wing_span>>1;
  360. CCshow (p, 2, "b2", DV(b2*100L));
  361.         b2VV = VONE*16*b2;
  362. CCshow (p, 0, "b2VV", (long)b2VV);
  363.         c2 = EP->wing_cord>>1;
  364. CCshow (p, 2, "c2", DV(c2*100L));
  365.         c2VV = VONE*16*c2;
  366. CCshow (p, 0, "c2VV", (long)c2VV);
  367.  
  368.         Crudder = -muldiv (EP->MaxRudder, EX->rudder, RAD2ANG(100));
  369.         Cailerons = stall ? 0 :
  370.             -muldiv (EP->MaxAilerons, EX->ailerons, RAD2ANG(100));
  371.         Celevators = stall ? 0 :
  372.             -muldiv (EP->MaxElevators, EX->elevators, RAD2ANG(100));
  373.         Cb = beta;
  374.         Cdx = p->da[X];
  375.         Cdy = p->da[Y];
  376.         Cdz = -p->da[Z];
  377. CCshow (p, 3, "Crdr", (long)fmul(Crudder, 1000));
  378. CCshow (p, 3, "Cail", (long)fmul(Cailerons, 1000));
  379. CCshow (p, 3, "Celev", (long)fmul(Celevators, 1000));
  380. CCshow (p, 3, "Cb", (long)fmul(Cb, 1000));
  381. CCshow (p, 3, "Cdx", (long)fmul(Cdx, 1000));
  382. CCshow (p, 3, "Cdy", (long)fmul(Cdy, 1000));
  383. CCshow (p, 3, "Cdz", (long)fmul(Cdz, 1000));
  384.  
  385.         t = fmul (Cl*CLf, vmagV>>1);
  386.         lift = muldiv (rVS, t, weight)*2;    /* cheat */
  387. CCshow (p, 0, "lift", (long)lift);
  388.         t = fmul (Cd, vmagV>>1);
  389.         drag = muldiv (rVS, t, weight);
  390. CCshow (p, 0, "drag", (long)drag);
  391.  
  392. /* This formula is faking behaviour at high aoas (usually when stall is
  393.  * disabled).
  394. */
  395. #if 0
  396.         drag += fmul (iabs(lift), FONE-ClStall);
  397. #endif
  398.         lift = fmul (lift, ClStall);
  399.  
  400. /* Fx,Fy,Fz: aerodynamic forces on the plane, operating on the aero. center.
  401. */
  402.         t =    fmul (vmagV, fmul (EP->Cydr, Crudder)) +
  403.             fmul (Cb, EP->Cybeta)*10
  404.             ;
  405.         Fx = muldiv (rVS, t, weight);
  406.         Fx += fmul (sb, ihypot2d (drag, lift));
  407. CCshow (p, 0, "Fx", (long)Fx);
  408.         Fy = fmul (cb, fmul (sa, lift) - fmul (ca, drag));
  409. CCshow (p, 0, "Fy", (long)Fy);
  410.         Fz = fmul (cb, fmul (ca, lift) + fmul (sa, drag));
  411. CCshow (p, 0, "Fz", (long)Fz);
  412.  
  413.         t = fmul (ClTail, vmagV>>1);
  414.         t = muldiv (rVS, t, weight);    /* Tail lift */
  415.         t = fmul (ca, t);
  416. CCshow (p, 0, "Tlift", (long)t);
  417.  
  418.         t =    fmul (c2VV, fmul (EP->Cmq, Cdx))*10    /* damping */
  419.               +    fmul (t, EP->Cmalpha)            /* stabilizers */
  420.               +    fmul (vmagV, fmul (EP->Cmde, Celevators));
  421.         t = muldiv (c2, t, VONE);
  422.         Mx = muldiv (rVS, t, VONE*VONE*VONE);
  423.  
  424. #define Clr    (Cl/4)
  425.  
  426.         t =    fmul (b2VV, fmul (EP->Clp, Cdy))    /* damping */
  427.               + fmul (b2VV, fmul (Clr, Cdz))
  428.               + fmul (EP->Clbeta, Cb)
  429.               + fmul (vmagV, fmul (EP->Clda, Cailerons));
  430.         t = muldiv (b2, t, VONE);
  431.         My = muldiv (rVS, t, VONE*VONE*VONE);
  432.  
  433.         t =    fmul (vmagV, fmul (EP->Cndr, Crudder))
  434.               + fmul (b2VV, fmul (EP->Cnp, Cdy))
  435.               + fmul (b2VV, fmul (EP->Cnr, Cdz))    /* damping */
  436.               + fmul (EP->Cnbeta, Cb)
  437.               +    fmul (vmagV, fmul (EP->Cnda, Cailerons));
  438.         t = muldiv (b2, t, VONE);
  439.         Mz = muldiv (rVS, t, VONE*VONE*VONE);
  440. CCshow (p, 0, "Mx", (long)Mx);
  441. CCshow (p, 0, "My", (long)My);
  442. CCshow (p, 0, "Mz", (long)Mz);
  443.  
  444.         AA[X] += Fx;
  445.         AA[Y] += Fy + push;
  446.         AA[Z] += Fz;
  447.  
  448. #if 0
  449. /* Examining the moment graphs at the end of chapter 13 of "The Design Of The
  450.  * Aeroplane" [Darrol Stinton, pub. Collins] one can roughly extract the
  451.  * following values for weights around 15000Kg (the very top end):
  452.  * Roll  moment (Ixx) = weight * 9.
  453.  * Pitch moment (Iyy) = weight * 7.
  454.  * yaw   moment (Izz) = weight * 17.
  455.  * Ixz is ignored.
  456.  * Note that our x,y,z are non-standard!
  457. */
  458.         Ixx = fmul (weight, FCON(0.72));
  459.         Iyy = fmul (weight, FCON(0.17));
  460.         Izz = fmul (weight, FCON(0.86));
  461. #else
  462.         Ixx = fmul (weight, EP->Iyy);
  463.         Iyy = fmul (weight, EP->Ixx);
  464.         Izz = fmul (weight, EP->Izz);
  465. #endif
  466. CCshow (p, 0, "Ixx", (long)Ixx);
  467. CCshow (p, 0, "Iyy", (long)Iyy);
  468. CCshow (p, 0, "Izz", (long)Izz);
  469.  
  470.         da[X] += muldiv (DEG(57)/VONE, Mx, Ixx);
  471.         da[Y] += muldiv (DEG(57)/VONE, My, Iyy);
  472.         da[Z] -= muldiv (DEG(57), Mz, Izz);    /* note *VONE */
  473.         if (stall) {
  474.             da[X] += Frand()%(VD90/16*2+1) - (VD90/16);
  475.             da[Y] += Frand()%(VD90/16*2+1) - (VD90/16);
  476.             da[Z] -= Frand()%(VD90/16*2+1) - (VD90/16);
  477.         }
  478.  
  479.         if (onground) {
  480.             t = DV(p->speed) - (EP->liftoff_speed>>1); /*nm->meter*/
  481.             if (t < 0 && p->a[X] <= EP->gpitch)
  482.                 t = 0;
  483.             if (da[X] > t)
  484.                 da[X] = t;
  485.             if (da[X] < 0 && p->a[X] <= EP->gpitch)
  486.                 da[X] = 0;
  487.  
  488.             if (p->a[Y])
  489.                 p->a[Y] -= TADJ(p->a[Y]); /* level wings */
  490.  
  491.             da[Z] += -TADJ(p->da[Z]); /* no rotation */
  492.  
  493.             p->da[Y] = 0;
  494.             da[Y] = 0;
  495.  
  496. /* Effective weight is reduced by lift (it is assumed that the lift is
  497.  * directed straight up which is reasonable when one is moving parallel to
  498.  * the ground). We must be going fast enough for the brakes to be unable to
  499.  * completely stop us within this interval.
  500. */
  501.             drag = EP->wheel_mu + muldiv (EP->brake_mu-EP->wheel_mu,
  502.                             EX->brake, 100);
  503.             drag = fmul (drag, lift - GACC);
  504.             if (EX->v[Y] < 0)
  505.                 drag = -drag;
  506.             AA[Y] += fmul (ca, drag);
  507.             AA[Z] += fmul (sa, drag);
  508.             EX->flags &= ~PF_STALL;
  509.             stall = 0;
  510.         } else {
  511.             EX->Gforce += AA[Z];
  512.             if (EX->Gforce > EP->max_lift || EX->Gforce < EP->min_lift)
  513.                 EX->flags |= PF_GLIMIT;
  514.             if (EX->equip & EQ_GEAR) {    /* gear shaking */
  515.                 t = ~1 & fmul (p->speed, FCON(0.02));
  516.                 t = t*t;
  517.                 AA[X] += TADJ(Frand()%(1+t) - (t>>1));
  518.                 AA[Y] -= TADJ(Frand()%(1+t));
  519.                 AA[Z] += TADJ(Frand()%(1+t) - (t>>1));
  520.             }
  521.         }
  522.  
  523.         EX->v[X] += TADJ(AA[X]);
  524.         EX->v[Y] += TADJ(AA[Y]);
  525.         EX->v[Z] += TADJ(AA[Z]);
  526.         p->speed = ihypot3d (EX->v);
  527.         VMmul (p->V, EX->v, p->T);
  528. #if 0
  529.         if (onground && p->V[Z] < 0) {
  530.             p->V[Z] = 0;
  531.             VxMmul (EX->v, p->V, p->T);
  532.         }
  533. #endif
  534.     }
  535.  
  536.     p->da[X] += TADJ (da[X])*VONE;    /* rotations */
  537.     p->da[Y] += TADJ (da[Y])*VONE;
  538.     p->da[Z] += TADJ (da[Z]);
  539.  
  540.     da[X] = TADJ(p->da[X])*VONE;
  541.     da[Y] = TADJ(p->da[Y])*VONE;
  542.     da[Z] = TADJ(p->da[Z])*VONE;
  543.     Myxz (p->T, da);
  544.  
  545.     if (!taxiing) {
  546.         Vcopy (AA, EX->v);
  547.         VxMmul (EX->v, AA, p->T);
  548.     }
  549.  
  550.     fMroty (p->T, p->siny, p->cosy);
  551.     fMrotx (p->T, p->sinx, p->cosx);
  552.     fMrotz (p->T, p->sinz, p->cosz);
  553.     Mangles (p, p->T, p->a, da[Y]);
  554.  
  555.     if (onground && p->a[X] < EP->gpitch) {
  556.         p->a[X] = EP->gpitch;
  557.         Mobj (p);
  558.     }
  559.  
  560.     if (taxiing) {
  561.         VMmul (p->V, EX->v, p->T);
  562.         p->V[Z] = 0;
  563.     }
  564.  
  565.     if (onground && p->V[Z] < 0) {
  566.         p->V[Z] = 0;
  567.         VxMmul (EX->v, p->V, p->T);
  568.         EX->v[X] = 0;    /* testing */
  569.         p->speed = ihypot3d (p->V);
  570.     }
  571.  
  572. #define MAX_SPEED    (1000*VONE)
  573.     if (p->speed > MAX_SPEED) {            /* temp */
  574.         t = fdiv (MAX_SPEED, p->speed);
  575.         EX->v[X] = fmul (t, EX->v[X]);
  576.         EX->v[Y] = fmul (t, EX->v[Y]);
  577.         EX->v[Z] = fmul (t, EX->v[Z]);
  578.         p->V[X]  = fmul (t, p->V[X]);
  579.         p->V[Y]  = fmul (t, p->V[Y]);
  580.         p->V[Z]  = fmul (t, p->V[Z]);
  581.         p->speed = ihypot3d (p->V);
  582.     }
  583. #undef MAX_SPEED
  584.  
  585. /* Mach number.
  586. */
  587.     EX->mach = muldiv (p->speed, 1000, sos);
  588.  
  589. /* pull up warning time.
  590. */
  591.     t = muldiv (4000, iabs(p->speed), 300*VONE);
  592.     t = muldiv (t, iabs(p->a[X]), D90);
  593.     if (t < 2000)
  594.         t = 2000;
  595.     EX->misc[8] = t;
  596. }
  597.  
  598. #undef CLf
  599. #undef Clr
  600.